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such  a  simple  function  of  z,  then  the  corresponding  eigenfunctions  can  be 
generated  from  the  above  set  by  procedures  analogous  to  the  Rayleigh- 
Schrodinger  perturbation  theory  used  in  quantum  mechanics.  The  avail¬ 
ability  of  such  eigenfunctions  enables  one  to  evaluate  quantitatively  the 
manner  in  which  geomagnetically  trapped  particles  are  redistributed  in 
<*0  and  lost  from  the  magnetosphere  as  the  phase-space  density  f  evolves 
in  time. 
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FIGURE 


1.  Domain  of  integration  for  evaluating  particle  content 

of  a  drift  shell,  as  in  (49 )- ( 5 1 ) .  27 
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INTRODUCTION 


It  is  a  well-known  result  in  magnetospheric  physics  that  the  phase- 
space  density?  (averaged  over  gyration,  bounce,  and  drift)  evolves  ac¬ 
cording  to  the  equation 


K 

at 


1 

x  T(y) 


9f  1 
Dxx  ax  J 


+  s 


(1) 


in  the  presence  of  (a)  pitch-angle  diffusion  at  fixed  particle  energy  E  and 
shell  parameter  L  and  (b)  a  distributed  source  S  (MacDonald  and  Walt, 

1961;  Haerendel,  1968;  Roberts,  1969;  Walt,  1970;  Lyons  et  al.  ,  1972; 
Schulz  and  Lanzerotti,  1974).  The  diffusion  coordinate  x  is  the  cosine 
of  the  equatorial  pitch  angle  in  this  formulation,  and  the  factor  T(y)  is 
well  approximated  (Davidson,  1976)  by  the  formula 

T(y)  *  T(0)  -  [T(0)  -  T(l)]y3/4,  (2) 

1/2 

where  T(0)  asl.  3801730,  T(l)  =  (ir/6)(2)  «  0,7404805,  and  y  =  sin 

2  1/2 

(1  -  x  )  .  Eigenfunction  solutions  of  (1)  have  been  obtained  by  MacDonald 

and  Walt  (1961)  and  by  Roberts  (1969)  for  particular  functional  forms  of  the 
bounce-averaged  diffusion  coefficient  D^  under  the  approximation  that  T(y) 
commutes  with  d/dx.  The  difficulty  with  such  an  approximation  is  that  it  is 
credible  only  for  x  «  1,  whereas  one  often  requires  solutions  that  are  valid 
over  the  entire  interval  Osx^sl. 


The  purpose  of  the  present  work  is  to  introduce  a  new  diffusion  co¬ 


ordinate  called  z  ,  in  terms  of  which  (1)  can  be  solved  without  further 
approximation  over  the  entire  range  of  x  for  selected  forms  of  Dzz  =s 
2 

[  xT(y)]  D  .  The  new  coordinate  is  defined  by  the  equation 
rx  r1 

z  =  Z(y)  =  /  x'T(y')dx'  =  /  y'  T(y')  dy' 

o  *y 

*  i.(l  -  y2)T(0)  -  -±.[t(0)  -  T(l)]  (1  -  y11/4)  (3) 

and  can  be  shown  (Schulz,  1974)  to  assume  the  end-point  values  Z(0) 

=  16/35  and  Z(l)  =  0  exactly.  The  approximation  for  Z(0)  extracted 
from  (3)  agrees  with  16/35  to  within  0.  1%  (Schulz,  1976).  It  follows  from 
(1)  and  (3)  that 


-!L  =  j_rD  -|L1  +  §  (4) 

at  az  L  m  dz  J 

with  Dzz  defined  as  above.  This  last  form  of  the  diffusion  equation  is 

canonical  in  the  sense  that  there  is  no  intervening  Jacobian  factor  that 

fails  to  commute  with  3/8z .  Thus,  if  Dzz  is  a  suitably  simple  function 

of  z,  then  one  can  specify  the  eigenfunctions  gn(z)  of  the  diffusion  operator 

in  closed  form  by  requiring  gn(*c)  to  vanish  for  some  positive  z£<  16/35. 

The  resulting  eigenfunctions  will  be  applicable  to  the  entire  physical  range 

(0  £  z  s  z  <  16/35)  of  the  new  canonical  diffusion  coordinate  z. 
c 
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If  D  is  not  precisely  of  such  a  form  that  yields  gn(z)  in  terms 

of  previously  studied  analytical  functions,  it  may  nevertheless  happen 

that  D  closely  resembles  in  form  a  diffusion  coefficient  D  for  which 
zz  -  zz 

the  eigenfunctions  gn(z)  are  known  exactly.  In  this  case  one  may  be  able 
to  use  the  gn(z)  as  a  basis  for  generating  the  gn(z)  by  means  of  pertur¬ 
bation  theory.  These  and  other  applications  of  the  coordinate  z  are  examined 
below. 


EXACT  EIGENFUNCTIONS 

The  functional  form  of  D  is  neither  well  known  nor  easily  (cf.  Lyons 

et  al.  ,  1972)  derived.  However,  the  construction  of  pitch-angle  eigenfunctions 

gn(z)  for  (4)  can  be  illustrated  quantitatively  if  one  arbitrarily  assumes  that 

D  =  (z/z _)a  D  *  where  <r  <  2  and  D  ,c  is  the  value  of  D  at  some  z  =  z 
zz  c  zz'  zz  zz  c 

<16/35  where  f  is  required  to  vanish.  There  is  a  precedent  for  this  type  of 

exercise  in  the  work  of  Roberts  (1969),  who  sought  solutions  of  (1)  for  D  = 

xx 

(  x/x  )^D^£  under  the  assumption  that  T(y)  would  commute  with  9/9x  .  The 
present  form  of  D  agrees  with  the  Roberts  (1969)  form  of  D  for  x^  «  1 
if  one  takes  £  =  2<r  -  2,  but  both  forms  are  equally  arbitrary. 

There  is  a  broader  purpose  behind  an  exercise  of  this  type.  By 

identifying  certain  functional  forms  of  D  for  which  the  eigenfunctions  of 

the  diffusion  operator  can  be  expressed  in  closed  form,  one  thereby  obtains 

a  complete  set  of  orthogonal  functions  on  the  interval  Os  z  s  z  .  The 

c 
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eigenfunctions  corresponding  to  a  somewhat  different  form  of 

can  be  expanded  in  terms  of  this  set,  and  the  expansion  coefficients  can 

be  determined  by  means  of  perturbation  theory. 

For  this  purpose,  it  proves  useful  (see  above)  to  adopt  the  notation 

D  for  any  special  form  of  D  that  leads  to  eigenfunctions  which  can  be 
z  z  z  z 

written  explicitly  in  closed  form.  It  is  logical  then  to  denote  the  eigen¬ 
values  of  the  diffusion  operator  as  X.  and  the  corresponding  eigenfunctions 
as  In(z)  when  has  such  a  special  form.  Thus,  in  the  present  context, 
one  is  considering  the  special  case  in  which  D  =  (z/z  )<r  D  where 

Z Z  C  ZZ  f 

<r<  2  and  D  "  is  the  value  of  D  at  z  =  z  . 

ZZ  ZZ  c 

Following  the  mathematical  methods  of  Roberts  (1969),  one  seeks 
solutions  of  the  eigenvalue  equation 


(d/dz)  [  Dzz  (dgn/dz)  ]  +  \ngn  =  0 


(5) 


in  the  form  g  (z)  «  z°w(PzY)  for  D  =  (z/z  J^D  *  .  Since  (5)  then 

11  ZZ  c  ZZ 

2 

reduces  to  Bessel' s  equation  for  Y  =  1  -  (cr/2),  a  =  (1  -<r)/2,  and  p  = 
z9 1  E>zzY^  follows  that  the  eigenfunctions  of  (5)  are  given  by 


g„U>  -  [(2-,)/zcJ1/2[VUvti)]-l 


v  vn  c 


(6) 
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where  K  (n  =  0,  1 ,  2,  .  .  .  )  is  the  nth  zero  of  the  ordinary  Bessel 
vn  — 

function  of  order  v  =  (cr  -  l)/(2  -  cr),  i.  e.  ,  where  ( K  )  =  0.  The 
corresponding  eigenvalue  is  given  by 


X.  =  (K  /  2  z)Z  (2-o  )Z  D  *, 
n  vn  c  zz 


and  the  normalization  of  (6)  has  been  chosen  so  that 


/  C  gn<z>gm<z>dz  =  6, 


1,  n  =  m 


0,  n  ^  m 


2  1/2 

Sinze  z  is  given  by  (3)  as  a  function  of  y  =  (l  -  x  )  ,  it  is  easy  enough 

to  plot  g  (z)  as  a  function  of  the  more  familiar  variable  x.  However,  no 
n 

immediate  purpose  would  be  served  by  such  a  plot,  since  the  form  of  Dzz 

leading  to  (6)  and  (7)  does  not  correspond  exactly  to  the  functional  form  of 

Dxx  postulated  by  Roberts  (1969)  except  in  the  limit  x  — 0;  nor  does  it 

correspond  to  the  functional  form  of  D  used  by  anyone  else.  Therefore, 

illustration  of  the  functional  form  of  g  (z)  is  deferred  for  now  and  is 

given  instead  in  the  accompanying  numerical-applications  paper  by 

Schulz  and  Boucher  (1981),  wherein  the  eigenfunctions  corresponding 

to  D  =  (x/x  )^D  *  are  estimated  successively  for  comparison  with  the 
xx  '  c  xx 

results  of  Roberts  (1969). 

The  form  D  =  (z/z  )°"d  *  considered  above  is  not  the  only 
zz  c  zz 

form  of  D  that  yields  exact  expressions  for  the  eigenfunctions  of 
z  z 

(5).  Another  form  of  D  that  leads  to  exact  eigenfunctions  is  the  form 


9 


D  =  (1  -  a^z^)  D  ®  ,  where  0  <  az  <1  and  D  ^  denotes  the  value 
zz  zz  c  zz 

of  Dzz  at  z  -  0.  This  form  of  D  converts  (5)  into  Legendre's  equation 
in  the  variable  t,  s  q-z>  Solutions  are  given  in  unnormalized  form  by  the 
expres  sion 


g  (z)  «  Q  (az  )  P  (az)  -  P  (az  )Q  (az), 

°V/  y  c  v  V  c  v 


(9) 


where  P^(f,)  and  Q^(£)  are  the  two  kinds  of  Legendre  function  of  degree  v 
and  order  p  =  0  (Stegun,  1966).  The  corresponding  eigenvalues  X.  are 


given  by 


X  =  aZ  v  (v  +  1)  D  °. 
v  zz 


(10) 


Acceptable  values  of  v{>0)  are  restricted  by  the  condition 


R  (az.  )  h  cos  (vn/Z)  P  (az  ) 
v  c  v  c 


-  (2/tt)  sin  ( vir/2)  Q^(  a z^)  =  0 ,  (11) 


which  assures  that  lim  [D  gv(z)]  =  0  as  z  -* 0.  The  significance  of 
this  latter  requirement  is  that  there  must  be  no  diffusion  current  across 
the  "boundary"  at  z  =  0,  since  particles  cannot  be  lost  by  having  their 
mirror  points  reach  the  equator.  The  eigenfunctions  gv(z)  and  gp(z) 
corresponding  to  distinct  eigenvalues  X^  and  X^  are  orthogonal  in  the 
sense  of  (8),  but  calculation  of  the  appropriate  normalization  constant 
for  (9)  in  closed  form  appears  to  be  intractable. 
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The  case  a  =  0  is  not  included  in  (9)- (11).  Eigenfunctions  in  this 


Jv(*vn>  ~  (ZhKVn)l'Z  C0SKn  '  v(Tr/2)  '  ' 

which  is  to  say  that 


(14) 


case  are  given  by 


1/2 


gn  (z)  =  (2/zc>  cos  [  (2n  +  l)(trz/2z^)  ], 


and  the  corresponding  eigenvalues  are  given  by 


X  =  (2n  +  l)2  (tt/2z  )2  D  0 
n  c  z  z 


(12) 


(13) 


for  n  =  0,1,2,...;  this  case  is  also  covered  by  (5)-  (8)  with  <r  =  0(v=-l/2), 

1/2 

since  one  knows  (e.  g.  ,  Antosiewicz,  1966)  that  J  j.gd*)  ~  (2/tt  4)  x 

cos  £.  Moreover,  one  finds  =  Dzz  by  definition  for  tr  =  0. 


The  distribution  of  eigenvalues  is  immediately  apparent  from  (13) 

for  the  case  in  which  D  is  independent  of  z.  However,  the  explicit  forms 

z  z 

of  k  in  (7)  and  of  v  in  (10)  are  available  only  in  the  asymptotic  (n  —  <o  ) 
limit  if  Dzz  varies  with  z.  Asymptotic  expansion  of  (*vn)  =  0  in  (6)  for 
large  argument  yields 


K  ~  nir  +  [(4  -  <r)/(2  -  «r)](ir/4). 
vn 


(15) 


The  asymptotic  expansion  of  If  given  by  (15)  yields  all  the  eigenvalues 
exactly  for  a  =  0  [see  (13)]  but  is  only  indicative  for  other  values  of  <r  <  2. 
On  the  other  hand,  the  asymptotic  expansion  of  R^faz^)  =  0  in  (11)  for 
large  v  yields 

R^az^l  ~  (2/ttv  sin0)^2cos|  [  v  +  (1/2)  ]  [8  -  (tt/Z) ]  }  =  0 ,  (16) 

where  8  s  cos~\arz  ).  This  means  that 
c 

v  ~  [  (2n  +  1)tt/(tt  -  28)  ]  -  (1/2)  (17) 

for  large  integers  n.  If  one  substitutes  (17)  in  (10)  and  takes  the  limit 
e-’-O,  the  eigenvalues  \vapproach  those  given  by  (13)  for  the  correspond¬ 
ing  values  of  n. 

PERTURBATION  THEORY 

It  would  be  fortuitous  if  the  form  of  D  in  a  realistic  situation 

zz 

corresponded  exactly  to  an  idealized  form  (D  )  known  to  yield  eigenfunc- 

tions  in  closed  form.  However,  it  is  not  unreasonable  to  expect  that  the 

physically  realistic  D  might  be  roughly  approximated  by  some  such  D 

zz  zz 
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In  this  case  one  might  be  able  to  use  the  eigenfunctions  g^(z)  that  correspond 

to  D  as  a  basis  (i.  e.  ,  as  a  complete  set  of  orthogonal  functions)  for  gen- 
z  z 

erating  by  means  of  perturbation  theory  the  eigenfunctions  g  (z)  that  cor¬ 
respond  to  D 
r  zz 

Let  the  linear  transformation  between  the  true  eigenfunctions  g^(z) 
and  the  basis  functions  g^(z)  be  specified  by 

co 

g_(z)  =  Y  g  (z)  U  ,  (18) 

n  / i °m  mn 

m  =  0 

where  the  expansion  coefficients  Umn  form  a  real  unitary  matrix,  i.  e.  ,  a 
matrix  such  that 


EU  U  =  6  .  (19) 

mp  mn  pn  v  ' 

m  =  0 

The  expansion  coefficients  Umn  are  otherwise  unknown  at  this  stage.  How¬ 
ever,  the  contention  that  gn(z)  is  an  eigenfunction  corresponding  to  &zz  must 
mean  that 

(d/dz)  [Dzz(dgn/dz)  ]  +  \ngn  =  0  (20) 

for  some  \  [compare  with  (5)  ].  By  invoking  the  orthogonality  property 
specified  by  (8),  one  thereby  derives  from  (18)  and  (20)  the  condition 

GO 

(21) 

pm  n  pm  mn 

m  =  0 
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on  the  expansion  coefficients  U  ,  where 

mn 


Apm  S  fC Ip(d/dz)[Dzz(dgm/dz)]  dz 


r  c  _  __ 

J  (dgp / dz)  Dzz  (dgm/dz)  dz. 


The  second  (i.  e.  ,  the  manifestly  symmetric)  integral  expression  for  A 
is  derived  from  the  first  through  integration  by  parts.  One  makes  use 


(22) 


pm 


here  of  the  fact  that  g  (z  )  =  0  and  the  requirement  (compare  above)  that 

P  c 

lim  [Dzz  gm  (z)  ]  =  0  as  z-*0.  This  latter  requirement  means  that  one  must 

select  D  so  that  lim  (D  /D  )  *  0  as  z -»  0,  i.e.  ,  so  that  the  g  (z)  do  not 
zz  zz  zz  em 

transport  particles  across  the  kinematical  "boundary"  at  z  -  0. 


It  follows  from  (21)  that  the  columns  of  the  unitary  matrix  U  are  the 

mn 

normalized  eigenvectors  of  the  real  symmetric  matrixA  ,  and  that  the  X. 

pm  n 

are  the  corresponding  eigenvalues.  Thus,  the  problem  of  identifying  the 
eigenfunctions  and  eigenvalues  of  (20)  has  been  reduced  to  the  problem  of  di¬ 
agonalizing  the  matrix  Apm  specified  by  (22).  One  observes  from  the  first 

integral  expression  for  A  that  if  D  =  D  then  A  =6  X.  .  This 

pm  zz  zz,  pm  pm  m 

follows  from  (5)  and  (8).  Consequently,  if  Dzz  is  a  reasonably  good  approxi¬ 
mation  of  D  then  the  off-diagonal  elements  of  A  will  be  "small"  in  the 
zz,  °  pm 


sense  required  by  perturbation  theory. 
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I 


The  usual  procedure  for  diagonalizing  a  matrix  such  as  Apm  is  to 
seek  (eigenvalue)  solutions  X.  of  the  characteristic  equation 


det(A  -6  X.)  =  0. 

pm  pm 


When  one  seeks  to  expand  the  above  determinant  by  minors,  it  becomes 
obvious  that  factors  lying  off  the  main  diagonal  can  affect  X  only  to  second 
or  higher  order.  Thus,  if  terms  of  second  and  higher  order  are  neglected, 
one  obtains 


A  «  A  =  X.  +  /*  (D  -D)  (dg  /dz)  dz 
n  nn  n  J  zz  zz  °n 

»  A 


for  the  eigenvalues  and 


uk„/u„„  *  Ak„/<A„„-Akk>'  k-” 


for  the  components  of  the  corresponding  eigenvectors.  One  can 
normalize  ( Z 5)  in  accordance  with  (19)  by  setting 


-  i‘  +  2>kn/<A„„-Akk>)2|-1/2 

k  *n 


although  it  is  apparent  from  (26)  that  U  «  1  except  for  corrections  of 

nn  r 

second  or  higher  order  in  "small"  quantities. 


The  foregoing  results  actually  represent  somewhat  of  an  im¬ 
provement  (cf.  Morse  and  Feshbach,  1953)  over  those  obtained  by 
the  usual  Rayleigh-Schrodinger  method  encountered  in  quantum  me- 
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chanics  (e.  g.  ,  Schiff,  1955).  The  advantage  of  the  present  method 
over  Rayleigh-Schrodinger  is  that  the  denominator  in  (25)  contains 
a  better  approximation  to  the  difference  between  the  true  eigenval¬ 
ues.  The  author  has  been  informed  by  Cornwall  (1977)  that  the 
present  perturbation  method  is  known  in  quantum  mechanics  as  the 
Wigner -Brillouin  method.  Expanding  (23)  to  second  order  in  off- 
diagonal  elements,  one  readily  obtains 


\ 


n 


A 


nn 


E 

k*n 


Akn Ank 
Akk  ^nn 


(27) 


as  an  improvement  on  the  second- order  Rayleigh-Schrodinger  result 
for  non-degenerate  states,  which  the  diffusion  eigenfunctions  clearly 
are.  Substitution  of  (27)  in  (21)  yields  a  general  equation  of  the  form 


A  + 

nn 


J*n 
+  A 


Aflj  Ajn 

Ajj  '  An» 


U. 

kn 

U 

nn 


,  +  y''  a,  .  (u.  /u  ) 

kn  /  7  kj  jn  nn 


i*k,n 


0 


(28) 


for  calculating  the  second-order  eigenvectors.  Since  A,  .  and  U.  are  both 

kj  jn 

"small"  quantities  for  k  *  j  *  n,  it  will  be  sufficient  (for  the  second-order 
accuracy  of  the  off-diagonal  elements  U^n)  to  estimate  the  ratios  U^/U^  in 
(28)  by  means  of  (25).  Similarly,  one  can  neglect  the  summation  that  appears 
in  the  square-bracketed  coefficient  of  U^/U^  in  (28)  without  sacrificing  the 
desired  order  of  accuracy.  Thus,  it  follows  from  (28)  that 


U 


kn 


U 


nn 


^nn  Akk 


Akn  + 


A,  .  A. 

kj  Jn 


A  -  A..  J 

j*k,n  nn  3J 


(29) 


16 


for  k*n.  The  diagonal  elements  Unn  are  to  be  determined  from  (19). 
A  first-order  expansion  of  (26)  assures  unit  normalization  to  second 
order  in  "small"  quantities.  However,  the  use  of  (29)  in  (19)  might 
be  preferable,  in  that  this  procedure  would  assure  unit  normalization 
of  each  perturbed  eigenvector  to  all  orders. 


WKB  APPROXIMATION 


An  alternative  construction  of  eigenfunctions  for  the  case  in 
which  D  varies  only  weakly  with  z  is  familiar  from  the  literature  of 
quantum  mechanics.  This  is  the  method  of  Wentzel  (1926),  Kramers 
(1926),  and  Brillouin  (1926).  The  WKB  approximation  is  motivated  by 
transforming  (20)  into  the  time-independent  Schrodinger  equation 

(d2gn/dC2)  +  k2gn  =  0.  (30) 


This  is  achieved  by  introducing  the  new  variable 

z 

1  -  J  <D«'D, 

0 

whereupon  one  obtains 


k2  =  \  (D  Y2D  • 
n  zz  zz 


(31) 


(32) 


Indeed,  equations  ( 30) - ( 32 )  are  valid  even  for  an  arbitrary  variation 
D zz  with  z.  However,  if  D  varies  only  weakly  with  z,  then  k  must 
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vary  only  weakly  with  £  .  In  this  case  one  obtains 


(33) 

as  the  WKB  solution  (e.  g.  ,  Merzbacher,  1961).  The  corresponding 
eigenvalues,  determined  by  requiring  that  gn(zc)  =  0,  are  given  by 


n 


a  (2n+  l)2  (tr/2)2Dz*  [  f 


(D  */D  )1/Z  dz 

zz  zz 


-2 


(34) 


The  above  results  for  g^(z)  and  reduce  (as  required)  to  (12)  and 
(13),  respectively,  in  case  D  is  a  true  constant  (altogether  indepen- 
dent  of  z).  They  provide  a  viable  alternative  to  perturbation  theory 
(based  on  the  above-described  case  cr  =  0)  if  D  varies  weakly  with  z. 

ZZ 


VARIATIONAL  PRINCIPLE 


A  further  point  deserves  consideration,  namely  that  the  diago- 
nalization  of  Apm  in  (21)  is  equivalent  to  the  implementation  of  a  vari¬ 
ational  principle  (Cornwall,  1977)  based  on  (22).  The  variational  prin¬ 
ciple  asserts  that 
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3 

lin  I 
J r\ 


X0  =  Min  /  (dg/dz)  D^z  dz 


subject  to  the  constraint  that 


zc 

J  [g(z)]‘ 


dz  =  1 , 


the  condition  that  lim  [D  g'(z)]  =  0  as  z-*-0,  and  the  boundary  condition 

zz 

that  g  <z^)  =  0. 

Proof:  Since  the  eigenfunctions  gn(z)  of  the  true  diffusion  operator 

(d/dz)  [D  (d/dz)]  form  a  complete  set  of  orthogonal  functions,  one  can 
z  z 

expand  any  continuous  g(z)  satisfying  (36)  so  that 


«<*>  =  £  v„Sn(a>- 
n  =  0 


where 


-  £  * 


i  ■***>*/* 


The  minimization  specified  by  (35)  is  implemented  by  varying  the  (real) 
expansion  coefficients  V  .  It  follows  from  (22)  ,  (37),  and  (38)  that 

z 

(dg/dz)2Dzz  dz 

CO 

=  \n  +  V*  (X  -  \n)  v2 .  (39) 

0  n  O  n 

n  =  l 

Since  Xn  >  Xq  for  n2  1,  i.e.,  since  one  of  the  eigenvalues  in  (20)  must 
be  the  smallest,  the  integral  that  appears  in  (35)  and  (39)  can  be 
minimized  only  by  setting  =  0  for  n  £  1.  Thus,  the  integral  in  (35) 
and  (39)  can  be  minimized  only  by  taking  g{z)  =  gg(z),  in  which  case  the 
integral  becomes  equal  to  Xq.  This  is  the  standard  proof  (e.g.,  Schiff,  1955) 
for  the  validity  of  a  variational  principle. 

The  usual  means  of  implementing  (35)  is  to  construct  a  trial 
function  g(z;c*m)  that  meets  the  required  constraints  and  depends  on  several 
adjustable  parameters  ar^.  The  integral  that  appears  in  (35)  is  then 
minimized  by  varying  the  adjustable  parameters.  It  is  not  always  practical 
to  obtain  eigenfunctions  higher  than  gg(z)  by  variational  means.  One 
theoretically  can  do  so,  as  in  quantum  mechanics  (e.g.  ,  Merzbacher,  1981), 
by  selecting  trial  functions  that  are  orthogonal  [in  the  sense  of  (8)]  to  the 
optimal  gQ(z)  obtained  in  the  manner  described  immediately  above.  Such  a 
procedure  is  often  too  cumbersome  for  practical  use,  especially  if  the 
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adjustable  parameters  °mn  appear  nonlinearly  in  the  specification  of 

the  trial  functions  e  (z;  a  ).  However,  a  special  case  of  the  variational 
6n  mn  r 

method  is  realized  if  one  specifies  each  g  (z;  a  )  as  a  linear  super- 

position  of  orthogonal  functions  gm(z),  as  in  (18).  In  this  case  the 

correspond  to  the  expansion  coefficients  Umn  in  (18),  and  minimization 

of  \  in  (35)  is  equivalent  to  diagonalization  of  the  matrix  as  in  (21). 

Of  course,  the  matrices  A  and  U  in  (21)  are  of  infinite  dimension. 

pm  mn 

This  precludes  their  numerical  evaluation  in  complete  form.  However, 
progressively  better  variational  approximations  of  the  eigenfunctions 
g^(z)  can  be  obtained  by  diagonalizing  progressively  larger  finite  sub¬ 
matrices  of  Apm,  i.e.  ,  by  truncating  the  summation  in  (21)  at  m  =  N-l 
for  progressively  larger  values  of  N. 
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STEADY  STATE 


If  the  phase -averaged  source  S  in  (4)  were  constant  in  time,  then 
the  solution  f  (z,  t)  would  approach  a  steady-state  solution  f^z) 
in  the  limit  t— ►  oo.  Following  Roberts  (1969),  one  can  obtain  this 
iw(z)  by  integrating  (4)  twice  with  respect  to  z  for  3f  /  9t  =  0. 

The  result 

zc  z' 

fj*>  =  /  (D,,,.)-1  /  S(z")dz"dz'  (40) 

z  *'o 

is  obtained  upon  application  of  the  relevant  boundary  conditions. 

Roberts  (1969)  has  noted  that  if  S  is  assumed  independent  of  x  (i.  e. , 
independent  of  z  in  the  present  context),  then  the  functional  form  of  f^ 
tends  to  resemble  that  of  the  lowest  eigenfunction  g^.  This  tendency  can 
be  made  understandable  by  expanding  S  in  (4)  as  a  linear  combination  of 
the  orthogonal  eigenfunctions  g^.  One  thereby  obtains 

®  fZc 

«.(■>  =  ^  Xn 1  gn^z)/  S(z’)  gn(z')  dz',  (41) 

where  S(z')  2  0  by  definition  (i.  e.,  S  is  a  source).  As  long  as  S(z') 
in  (41)  is  constant  (or  at  least  relatively  structureless)  over  the  interval 
0  <  z  <  z  ,  it  is  likely  that  the  moment  of  S  with  respect  to  the 
positive -definite  eigenfunction  gg(z')  will  exceed  those  with  respect  to 
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the  oscillatory  (not  positive-definite)  eigenfunctions  gn(z  )  for  n  >  1. 

This  property  and  the  weighting  by  X.  ^  in  (4  1  )  would  account  for 
the  tendency,  noted  by  Roberts  (1969),  for  f^(z)  to  resemble  gQ(z)  in 
functional  form.  Of  course,  it  follows  from  the  completeness  of  the  gn(z) 
as  an  orthonormal  set  of  basis  functions  that  f^  will  coincide  exactly 
with  in  functional  form  only  if  S  itself  is  directly  proportional  to  gQ(z). 


OMNIDIRECTIONAL  FLUX 


It  would  be  appropriate  to  relate  the  formal  results  obtained  above 
to  physically  observable  quantities.  Consider  an  off- equatorial  point,  i.  e.  , 
one  at  which  the  local  magnetic -field  intensity  B  exceeds  the  equatorial 
value  Bq  on  a  field  line  identified  by  the  dimensionless  label  L.  It  is  well 
known  (e.  g.  ,  Schulz  and  Lanzerotti,  1974)  that  the  unidirectional  flux  of 
particles  (per  unit  energy  and  solid  angle  at  local  pitch  angle  a  )  is  equal 
to  p^f,  with  f  evaluated  at 


y  =  (Bq/B)1  sin  a. 


In  order  to  specify  the  omnidirectional  flux  (per  unit  energy)  at  this 
point  in  space,  one  must  integrate  p^f  over  the  unit  sphere  in  momentum 


(p)  space.  Thus,  it  follows  from  (42)  that 
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4tt 


=  4rrp 


i 

•T  r\ 


COS  O 


f  d(coso) 


B0/S 


/"  O' 

=  2irp2y  (B/B0)[l  -  y2(B/Bo)]'1/2  f  d(y2) 
Bq/B 

=  -2*pZJ  [1  -  y2(B/BQ)]1/2  Of/9z)  T(y)  d(y2) 


(43) 


at  B/Bg  2  1.  The  final  line  of  (43)  results  from  integration  by  parts 
and  serves  to  simplify  the  required  numerical  quadrature. 

The  phase- space  density  f  that  appears  in  (43)  can  be  expanded 
as  a  weighted  series  of  eigenfunctions  of  the  diffusion  operator  (Roberts, 

1969): 


f  -  fjz)  +  ^  An(E,L;t)  gn(z).  (44) 

n 

The  omnidirectional  flux  described  by  (43)  can  thus  be  written  in  the 
form 


J4„  '  2’P2  [<5/*cD*«>G.(B/B0> 

+  An(E.L;t)  Gn(B/BQ)],  (45) 

n 

where 

VB 

Gn<B/B0)  -  -f  [1  -  y2(B/B0)]1/2  gn(z)  T(y)  d(y2)  (46) 

Jyc 
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and 


B0/B 

GJB/Bq)  =  -(x^D^/S)  J  [1  -  y2(B/B0)]1/2  f>)  T(y)  d(y2). 

yc 

(47) 


Another  physical  quantity  of  interest  is  the  particle  content  C  (per 

unit  L  and  energy)  of  a  magnetic  drift  shell.  To  obtain  this,  one  must 

integrate  (1/v)  J  over  the  drift-shell  volume,  where  v  is  the  speed 

of  the  particle.  The  cross-sectional  area  of  an  infinitesimal  drift  shell 

2 

of  "width"  dL  is  2tr  La  (B^/B)  dL,  where  a  is  the  radius  of  the  earth. 
Therefore,  the  volume  per  unit  "width"  is  given  by 

TT 

(BQ/B)  (ds/d0)  de 

/r+1 

=  2irL2a3  /  sin60  d(cos  0)  =  4w(l  6/35)L2  a3,  (48) 

*'-1 


where  0  is  the  magnetic  colatitude  and  s  is  the  coordinate  that  measures 
arc  length  along  the  dipolar  field  line.  It  follows  from  the  above  consider¬ 
ations  and  from  (43)  that 


C  = 


4tt2  La2  (p2/v) 


Brt/B 


J  (ds/d0)  J  i  sec  a  d  (y2)  d0. 


(49) 
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2 

The  integral  over  y  in  (49)  contains  contributions  from  all  particle 
trajectories  that  mirror  at  a  higher  latitude  than  the  point  identified  by 
the  local  B/Bq.  The  domain  of  integration  is  illustrated  in  Figure  1. 

If  the  order  of  integration  is  reversed,  then  one  obtains 

(1/v  cos  a  )  ds  d(y2)  (50) 

upon  recalling  that  f  satisfies  Liouville's  theorem,  i.  e.  ,  that  f  depends 
on  y  but  not  on  s  .  The  upper  limit  in  (50)  represents  the  arc  length 
from  the  equator  (s  =  0)  to  the  mirror  point  of  a  particle  whose  equatorial 
pitch  angle  is  sin'^y. 

The  integral  of  (1/v  cos  o)  with  respect  to  s  in  (50)  represents  half 
the  bounce  period  of  the  particle,  i.e.,  is  equal  to  (ZLa/v)  T(y),  where 
T(y)  is  the  dimensionless  function  that  appears  in  (1)  and  (2).  Thus,  it 
follows  from  (50)  that 

f  y  T(y)  dy 

T  dz.  (51) 


C  =  16  ir2  L2  a3  ( p2 / v )  J 


1 


=  16  tt2  L2  a3  (p2/v)  f 

J  r\ 


C  =  4  it  La 


L  S 


+  s 


m 
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Figure  1.  Domain  of  integration  for  evaluating  particle  content 


This  relationship  brings  out  the  true  physical  significance  of  the  new 
diffusion  coordinate  z,  namely  that  z  is  a  direct  measure  of  phase-space 
volume  at  fixed  E  and  L.  Such  a  finding  comes  as  no  surprise;  if  it 
were  not  so,  then  z  would  not  be  a  canonical  coordinate  in  the  sense  of 
(4)  and  some  Jacobian  factor  would  intervene  there. 


FORMAL  CONSIDERATIONS 

The  canonical  coordinate  z  introduced  above  is  a  variable 
corresponding  to  the  equatorial  pitch  angle  Oq.  More  generally,  one 
may  wish  to  identify  canonical  coordinates  corresponding  (respectively) 
to  kinetic  energy  E  and  shell  parameter  L,  so  that  (in  the  presence  of 
energy  transport  and  radial  diffusion)  the  Fokker-Planck  equation  can 
be  written  in  the  canonical  form 


11 

9t 


+ 


3  3 

V  V  _3L  D 

La  dQ.  ij  ao. 

i=i  j=i  1  J 


+  s, 

(52) 


where  <  dQ^/dt  >  and  D„  are  the  transport  coefficients.  This  will  be 
the  case  if  the  new  coordinates  are  related  by  a  canonical  transformation 
(Goldstein,  1950)  to  the  three  adiabatic  invariants 
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I 

i 


M  =  (p2/2mQ  Bq)  (1  -  x2) 

(53) 

i 

J  -  2  Lap  Y(y) 

(54) 

j 

$  =  (2tt  a2/L)  (Bq  L3  )  sgn  q 

(55) 

of  charged-particle  motion  in  a  dipole  field,  i.  e.  ,  if  the  Jacobian 
9(M,  J,  <t>)  /  9(0^  Q^,  Qj)  of  the  transformation  from  (M,  J,  <t>)  to 
(Qj.  Q2,  Q3)  is  a  constant  (Haerendel,  1968).  The  particle  described 
by  (53)  —  (55)  has  charge  q  and  rest  mass  mQ.  Its  scalar  momentum 
is  p  =  (  y2  -  lJ^^mgC,  where  y  =  1  +  (E/mQc2)  and  c  is 
the  speed  of  light.  The  particle  described  by  (53)-(55)  executes  a 
drift  shell,  bearing  the  dimensionless  label  L,  on  which  the  equatorial 
magnetic  field  is  BQ  (proportional  to  L"  ).  The  particle  has  an  equatorial 
pitch  angle  aQ  =  sin  *y  =  cos_1x,  and  the  earth  has  a  radius  denoted  a. 
The  function  Y(y)  in  (54)  is  given  (Schulz,  1971;  Davidson,  1976)  by 

1 

Y(y)  *  2y  J  (y')~Z  T(y')  dy' 
y 

«  2T(0)  +  [6T(0)  -  8T(l)]y  -  8  [T(0)  -  T(l)]y3/4,  (56) 

where  T(y),  the  function  specified  by  (2),  is  equal  to  (p/4 Lay  mQ)  times 
the  particle's  bounce  period. 


i 

i 

f 

i' 
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The  transformation  from  (M,  J,  $)  to  (E,  x,  L)  has  a  Jacobian 


that  is  given  by 


3 (M,  J,  #)/  3(E,  x,  L)  = 


8iTa^  m. 


c  L2  v  ( 


V  -  1 


,1/2 


x  T (y)  sgn  q. 


(57) 


This  is  therefore  not  a  canonical  transformation,  since  its  Jacobian  depends 

on  all  three  of  the  new  kinematic  variables:  on  E  through  the  factor 
2  1/2 

y(  y  -  1)  ,  on  x  through  the  factor  x  T(y),  and  on  L  through  the 

2 

factor  L  .  Making  use  of  this  factorization,  however,  one  clearly  can 
construct  a  set  of  new  variables 


such  that  the  Jacobian  9(M,  J,  *)/  3(Qj,  Q2,  Q3)  is  indeed  a  constant. 

2  3/2  3  3 

Thus,  the  coordinates  (y  -  1)  as  (p/mQc)  ,  z  s  Z(y),  and  L  are 

canonical  in  the  present  sense  and  are  therefore  eligible  for  use  in  (52), 

of  which  (4)  is  a  special  case.  (Magnetospheric  electric  fields  related 

to  convection  and  corotation  are  implicitly  neglected  in  the  present 

work,  as  are  day-night  asymmetries  in  the  magnetic  field.  This 

simplifying  assumption  is  important  for  the  validity  of  the  above  Q^, 
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as  new  variables,  since  it  would  be  a  mistake  to  adopt  new  varia¬ 
bles  that  fail  to  remain  constant  around  an  adiabatic  drift  shell.  ) 

The  derivation  of  (4)  from  (52)  is  achieved  by  letting  all  of  the 

transport  coefficients  except  D  (s  D,,)  vanish.  In  a  description 

z  z  c*  c* 

based  on  the  uncanonical  variables  (E,  x,  L),  this  condition  would  be 
expressed  by  letting  all  of  the  transport  coefficients  except  vanish. 

In  this  description,  however,  one  must  insert  in  (52)  the  Jacobian  of  the 
transformation  from  the  canonical  action  variables  (M,  J,  9)  to  the 
uncanonical  variables  (E,  x,  L)  selected  for  their  conceptual  convenience 
(Haerendel,  1968;  Schulz  and  Lanzerotti,  1974).  One  thereby  obtains 


(61) 


where  =  E,  =  x,  =  L,  and  G  is  the  Jacobian  given  by  (57). 

One  obtains  (1)  from  (61)  by  letting  all  of  the  transport  coefficients 
except  vanish.  This  is  the  derivation  described  by  Haerendel  (1968). 

It  is  contingent  upon  the  fact  that  all  three  of  the  adiabatic  invariants 
(M,  J,  9)  are  canonical  action  variables  (e.  g.  ,  Schulz  and  Lanzerotti, 
1974).  Their  conjugate  phases  (angle  variables  ?.)  are  cyclic  coordinates 
for  the  unperturbed  Hamiltonian  of  charged -particle  motion  in  the  adiabatic 
guiding-center  approximation. 
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An  earlier  derivation  of  (1)  by  MacDonald  and  Walt  (1961)  had 
been  based  on  arguments  of  the  type  advanced  in  (48)-(51),  namely 
that  xT(y)  dx  is  a  direct  measure  of  phase-space  volume  at  fixed  E 
and  L.  Neither  they  nor  Roberts  (1969),  however,  chose  to  exploit 
the  coordinate  z  as  a  natural  variable  for  the  construction  of  eigen¬ 
functions.  They  chose  instead  to  let  T(y)  commute  with  9/9x  in  (1). 


SUMMARY 

The  major  point  of  the  present  work  has  been  to  introduce  the  new 

canonical  variable  z,  as  defined  by  (3),  in  order  to  simplify  the 

description  of  pitch-angle  diffusion  in  a  dipolar  magnetic  field.  Various 

applications  seem  to  follow  quite  naturally.  For  example,  one  can 

calculate  eigenfunctions  of  the  diffusion  operator  (9/9z)  [D  (  9/9z)] 

zz 

by  means  of  a  quantum-mechanical  perturbation  theory,  if  not  in  closed 
form.  The  availability  of  such  eigenfunctions  enables  one  to  calculate 
properly  the  temporal  evolution  of  the  phase-space  density  f  from  an 
arbitrary  initial  configuration  toward  an  asymptotic  steady  state.  It 
would  surely  be  possible  to  identify  further  applications  of  the  canonical 
diffusion  coordinate  z.  Those  noted  above  should  suffice  to  establish 
the  usefulness  of  the  scheme.  Numerical  results  illustrating  implementa¬ 
tion  of  the  methods  described  above  are  given  in  an  accompanying  paper 
(Schulz  and  Boucher,  1981).  Other  applications  are  left  (for  now,  at  least) 
to  the  imagination  of  the  reader. 
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LABORATORY  OPERATIONS 

The  Laboratory  Operations  of  The  Aerospace  Corporation  is  conducting 
experimental  and  theoretical  investigations  necessary  for  the  evaluation  and 
application  of  scientific  advances  to  new  military  concepts  and  systems.  Ver¬ 
satility  and  flexibility  have  been  developed  to  a  high  degree  by  the  laboratory 
personnel  in  dealing  with  the  many  problems  encountered  in  the  nation's  rapidly 
developing  space  and  missile  systems.  Expertise  in  the  latest  scientific  devel¬ 
opments  is  vital  to  the  accomplishment  of  tasks  related  to  these  problems.  The 
laboratories  that  contribute  to  this  research  are: 

Aerophysics  Laboratory:  Launch  and  reentry  aerodynamics,  heat  trans¬ 
fer,  reentry  physics,  chemical  kinetics,  structural  mechanics,  flight  dynamics, 
atmospheric  pollution,  and  high-power  gas  lasers. 

Chemistry  and  Physics  Laboratory:  Atmospheric  reactions  and  atmos- 
pheric  optics,  chemical  reactions  in  polluted  atmospheres,  chemical  reactions 
of  excited  species  in  rocket  plumes,  chemical  thermodynamics,  plasma  and 
laser-induced  reactions,  laser  chemistry,  propulsion  chemistry,  space  vacuum 
and  radiation  effects  on  materials,  lubrication  and  surface  phenomena,  photo¬ 
sensitive  materials  and  sensors,  high  precision  laser  ranging,  and  the  appli¬ 
cation  of  physics  and  chemistry  to  problems  of  law  enforcement  and  biomedicine. 

Electronics  Research  Laboratory:  Electromagnetic  theory,  devices,  and 
propagation  phenomena,  including  plasma  electromagnetics;  quantum  electronics, 
lasers,  and  electro- optics;  communication  sciences,  applied  electronics,  semi¬ 
conducting,  superconducting,  and  crystal  device  physics,  optical  and  acoustical 
imaging;  atmospheric  pollution;  millimeter  wave  and  far- infrared  technology. 

Materials  Sciences  Laboratory;  Development  of  new  materials;  metal 
matrix  composites  and  new  forms  of  carbon:  test  and  evaluation  of  graphite 
and  ceramics  in  reentry;  spacecraft  materials  and  electronic  components  in 
nuclear  weapons  environment;  application  of  fracture  mechanics  to  stress  cor¬ 
rosion  and  fatigue -induced  fractures  in  structural  metals. 

Space  Sciences  Laboratory:  Atmospheric  and  ionospheric  physics,  radia¬ 
tion  from  the  atmosphere,  density  and  composition  of  the  atmosphere,  aurorae 
and  airglow;  magneto  spheric  physics,  cosmic  rays,  generation  and  propagation 
of  plasma  waves  in  the  magnetosphere;  solar  physics,  studies  of  solar  magnetic 
fields;  space  astronomy,  x-ray  astronomy;  the  effects  of  nuclear  explosions, 
magnetic  storms,  and  solar  activity  on  the  earth's  atmosphere,  ionosphere,  and 
magnetosphere;  the  effects  of  optical,  electromagnetic,  and  particulate  radia¬ 
tions  in  space  on  space  systems. 
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